/* Appliation of velocity boundary condition. */
# include "scene.h"
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

# ifdef _OMP_VEL_BC_
# include <omp.h>
# endif

extern double inlet_velocity(double r, double t, short index);

extern void vel_bc(struct grid *g, short tflag, unsigned short iRK)
  {
   short nr, nth, nz, i, _i, j, _j, k;
   double ***u, ***v, ***w, *r, *z, *th, t, dv, rdu;

   if( (iRK > 0 && tflag != 0) || (iRK == 0 && tflag != 1) || (iRK > 2) )
       {
        fprintf(g->lptr, "\nERROR: In 'vel_bc()': (tflag, iRK) = (%d, %d).\nThis function is designed to handle only following cases:\n\t[tflag = 0, iRK = 1, 2] and [tflag = 1, iRK = 0]\n", tflag, iRK);
        exit(-1);
       }

   th = g->th; nth = g->nth;
   t = g->t + tflag*g->dt;


/* ATTENTION:: u[NrG-2][j][k], v[i][NthG-2][k], w[i][j][NzG-2] is required for two things:
      1. Evaluation of Neumann Pressure-BC, and
      2. Evaluation of RHS of the PPE.
   Therefore, these must be assigned everywhere, be it a solid boundary or otherwise. */

/* --------------- Applying BC for the SECTION - 0 ------------------ */
   u = g->u[0]; v = g->v[0]; w = g->w[0]; nr = g->nr[0]; nz = g->nz[0]; r = g->r[0]; z = g->z[0];

# ifdef _OMP_VEL_BC_
  #pragma omp parallel private(i,j,k,_i,_j,dv,rdu)
# endif
    {
  /* At the top-interface open to the section-1. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG+nz] = g->u[1][i][j][NzG];
          v[i][j][NzG+nz] = g->v[1][i][j][NzG];
          w[i][j][NzG+nz] = g->w[1][i][j][NzG];
         }

  /* At the side-walls, excluding the cell touching the top-section 1. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(k = NzG; k <= NzG+nz-2; k++) for(j = NthG; j < NthG+nth; j++)
         {
          cm = &g->icmx_rcs[0];

          q1 = q2 = 0.0;
          for(indx = 0; indx <= cm->iend[NrG+nr-cm->respadi]-cm->istart[NrG+nr-cm->respadi]; indx++)
             if(cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx == NrG+nr) coeff = cm->mtrx[NrG+nr-cm->respadi][indx];
             else{
                  q1 += cm->mtrx[NrG+nr-cm->respadi][indx]*v[cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx][j][k];
                  q2 += cm->mtrx[NrG+nr-cm->respadi][indx]*w[cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx][j][k];
                 }

          v[NrG+nr][j][k] = -q1/coeff;
          w[NrG+nr][j][k] = -q2/coeff;

       /* The normal component "u" would be determined from the divergence-free condition. */
          

          u[NrG+nr][j][k] = (r[NrG+nr-1]/r[NrG+nr+1])*u[NrG+nr-2][j][k];
//           v[NrG+nr][j][k] = -v[NrG+nr-1][j][k];
//           w[NrG+nr][j][k] = -w[NrG+nr-1][j][k];
         }


  /* At the side-wall for the ghost cell touching the top-section 1 */
     i = NrG+nr; k = NzG+nz-1;
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(j = NthG; j < NthG+nth; j++)
         {
       /* Because of the periodic BC in theta-direction. */
          dv = ((j == NthG) ? v[i-1][NthG][k]-v[i-1][NthG+nth-1][k] : v[i-1][j][k]-v[i-1][j-1][k]);
          u[i][j][k] = (1 - r[i]/r[i+1])*(-0.5*(r[i]+r[i+1])*w[i-1][j][k-1]/(z[k+1]-z[k]) + dv/(th[j+1]-th[j]));
          v[i][j][k] = -v[i-1][j][k];
          w[i][j][k] = 0;
         }

  /* At the INLET */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(i = NrG; i < NrG + g->nr_nozzle; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG-1] = -u[i][j][NzG];
          v[i][j][NzG-1] = -v[i][j][NzG];

          switch(iRK)
               {
                case 0: w[i][j][NzG-1] = inlet_velocity(0.5*(r[i]+r[i+1]), t, 0);
                        break;

                case 1: w[i][j][NzG-1] = inlet_velocity(0.5*(r[i]+r[i+1]), t, 0)\
                             + c2*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 1);
                        break;

                case 2: w[i][j][NzG-1] = inlet_velocity(0.5*(r[i]+r[i+1]), t, 0)\
                             + (a31+a32)*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                             + a32*c2*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2);
                        break;
               }

          w[i][j][NzG-2] = w[i][j][NzG];
         }

  /* At the bottom SOLID-BASE beyond the nozzle. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(i = NrG+g->nr_nozzle; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG-1] = -u[i][j][NzG];
          v[i][j][NzG-1] = -v[i][j][NzG];
          w[i][j][NzG-1] = 0.0;
          w[i][j][NzG-2] = w[i][j][NzG];
         }

  /* At the corner-cell after the SOLID-BASE. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic,3)
#    endif
     for(j = NthG; j < NthG+nth; j++) w[NrG+nr][j][NzG-1] = 0.0;

  /* For the NEGATIVE radius: applying symmetric B.C. across the center. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic,3)
#    endif
/*** The i-variation should always be kept in the end. ***/
     for(k = 0; k < 2*NzG+nz; k++) for(j = NthG; j < NthG+nth; j++) for(i = NrG; i >= 1; i--)
         {
     /*** Order of this if-else should not be reversed. ***/
          u[NrG-i][j][k] = (i > 1 ? -u[NrG+i-2][j][k]:0.5*(u[NrG-2][j][k]+u[NrG][j][k]));
          v[NrG-i][j][k] = -v[NrG+i-1][j][k];
          w[NrG-i][j][k] = w[NrG+i-1][j][k];
         }


  /* For ghost cells in theta-direction. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic, 2)
#    endif
  /* It needs the other boundary conditions applied first. */
     for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
      /* Applying Periodic B.C. */
         u[i][NthG-1-j][k] = u[i][NthG+nth-1-j][k];
         u[i][NthG+nth+j][k] = u[i][NthG+j][k];

         v[i][NthG-1-j][k] = v[i][NthG+nth-1-j][k];
         v[i][NthG+nth+j][k] = v[i][NthG+j][k];

         w[i][NthG-1-j][k] = w[i][NthG+nth-1-j][k];
         w[i][NthG+nth+j][k] = w[i][NthG+j][k];
        }


  /* --------------Applying BC for the top SECTION-1---------------- */
#    ifdef _OMP_VEL_BC_
     #pragma omp single
#    endif
       {
        u = g->u[1]; v = g->v[1]; w = g->w[1]; nr = g->nr[1]; nz = g->nz[1]; r = g->r[1]; z = g->z[1];
       }

  /* At the TOP-SURFACE */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG+nz] = -u[i][j][NzG+nz-1];
          v[i][j][NzG+nz] = -v[i][j][NzG+nz-1];
          w[i][j][NzG+nz] = w[i][j][NzG+nz-2];
         }

  /* At the SIDE-WALL */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < NthG+nth; j++)
         {
          u[NrG+nr][j][k] = (r[NrG+nr-1]/r[NrG+nr+1])*u[NrG+nr-2][j][k];
          v[NrG+nr][j][k] = -v[NrG+nr-1][j][k];
          w[NrG+nr][j][k] = -w[NrG+nr-1][j][k];
         }

  /* At the interface connected to the extended outlet-section 2. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(_i = 0; _i < g->nr_gap; _i++) for(i = _i+NrG+nr-g->nr_gap, j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG-1] = g->u[2][_i+NrG][j][NzG+g->nz[2]-1];
          v[i][j][NzG-1] = g->v[2][_i+NrG][j][NzG+g->nz[2]-1];
          w[i][j][NzG-1] = g->w[2][_i+NrG][j][NzG+g->nz[2]-1];
          w[i][j][NzG-2] = g->w[2][_i+NrG][j][NzG+g->nz[2]-2];
         }

  /* At the corner-cell after the interface connected to the extended outlet-section 2. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(j = NthG; j < NthG+nth; j++) w[NrG+nr][j][NzG-1] = -g->w[2][NrG+g->nr[2]-1][j][NzG+g->nz[2]-1];

  /* At the solid-surface between the connecting-interface to section-1 and the extended-outlet
     section-3 (excluding the end-cells touching the section-0 and section-2). */

#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(i = NrG+g->nr[0]+1; i <= NrG+nr-g->nr_gap-2; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG-1] = -u[i][j][NzG];
          v[i][j][NzG-1] = -v[i][j][NzG];
          w[i][j][NzG-1] = 0;
          w[i][j][NzG-2] = w[i][j][NzG];
         }

  /* At the cell touching the extended section-2. */
     i = NrG+nr-g->nr_gap-1; k = NzG-1;
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][k] = 0;
          v[i][j][k] = -v[i][j][k+1];
          w[i][j][k] = 0;

       /* Making use of the periodic BC in the theta-direction. */
          dv = (j == NthG ? v[i][NthG][k+1]-v[i][NthG+nth-1][k+1]:v[i][j][k+1]-v[i][j-1][k+1]);

          w[i][j][k-1] = 2*(z[k+1]-z[k])*(r[i]*u[i-1][j][k+1]/(r[i+1]-r[i]) - dv/(th[j+1]-th[j]))/(r[i+1]+r[i]);
         }

  /* At the cell touching the section-0. */
     i = NrG+g->nr[0]; k = NzG-1;
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][k] = -u[i][j][k+1];
          v[i][j][k] = -v[i][j][k+1];
          w[i][j][k] = 0;

       /* Making use of the periodic BC in the theta-direction. */
          dv = (j == NthG ? v[i][NthG][k+1] - v[i][NthG+nth-1][k+1] : v[i][j][k+1]-v[i][j-1][k+1]);

          w[i][j][k-1] = -2*(z[k+1]-z[k])*(r[i+1]*u[i][j][k+1]/(r[i+1]-r[i]) + dv/(th[j+1]-th[j]))/(r[i+1]+r[i]);
         }


  /* At the interface connected to the lower-section 0. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic,3)
#    endif
     for(i = NrG; i < NrG+g->nr[0]; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG-1] = g->u[0][i][j][NzG+g->nz[0]-1];
          v[i][j][NzG-1] = g->v[0][i][j][NzG+g->nz[0]-1];
          w[i][j][NzG-1] = g->w[0][i][j][NzG+g->nz[0]-1];
          w[i][j][NzG-2] = g->w[0][i][j][NzG+g->nz[0]-2];
         }

  /* For the NEGATIVE radius: applying symmetric B.C. across the center. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic,3)
#    endif
/*** The i-variation should always be kept in the end. ***/
     for(k = 0; k < 2*NzG+nz; k++) for(j = NthG; j < NthG+nth; j++) for(i = NrG-1; i >= 0; i--)
         {
     /*** Order of this if-else should not be reversed. ***/
          u[NrG-i-1][j][k] = (i > 0 ? -u[NrG+i-1][j][k]:0.5*(u[NrG-2][j][k] + u[NrG][j][k]));
          v[NrG-i-1][j][k] = -v[NrG+i][j][k];
          w[NrG-i-1][j][k] = w[NrG+i][j][k];
         }

  /* For ghost cells in theta-direction. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic,2)
#    endif
  /* It needs the other boundary conditions applied first. */
     for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
      /* Applying Periodic B.C. */
         u[i][NthG-1-j][k] = u[i][NthG+nth-1-j][k];
         u[i][NthG+nth+j][k] = u[i][NthG+j][k];

         v[i][NthG-1-j][k] = v[i][NthG+nth-1-j][k];
         v[i][NthG+nth+j][k] = v[i][NthG+j][k];

         w[i][NthG-1-j][k] = w[i][NthG+nth-1-j][k];
         w[i][NthG+nth+j][k] = w[i][NthG+j][k];
        }


 /* ------------Applying BC for the extende-outlet SECTION-2----------------- */
#    ifdef _OMP_VEL_BC_
     #pragma omp single
#    endif
        {
         u = g->u[2]; v = g->v[2]; w = g->w[2]; nr = g->nr[2]; nz = g->nz[2]; r = g->r[2]; z = g->z[2];
        }

  /* At the inner side-wall, excluding the top-cell touching the section-1. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(k = NzG; k <= NzG+nz-2; k++) for(j = NthG; j < NthG+nth; j++)
         {
          u[NrG-2][j][k] = (r[NrG+1]/r[NrG-1])*u[NrG][j][k];
          u[NrG-1][j][k] = 0;
          v[NrG-1][j][k] = -v[NrG][j][k];
          w[NrG-1][j][k] = -w[NrG][j][k];
         }

  /* At the inner side-wall for the top-cell touching the section-1. */
     i = NrG-1; k = NzG+nz-1;
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][k] = 0;

       /* Making use of the periodic BC in the theta-direction. */
          dv = (j == NthG ? v[i+1][NthG][k]-v[i+1][NthG+nth-1][k] : v[i+1][j][k]-v[i+1][j-1][k]);

          u[i-1][j][k] = (1 - r[i+1]/r[i])*(dv/(th[j+1]-th[j]) - 0.5*(r[i+1]+r[i])*w[i+1][j][k-1]/(z[k+1]-z[k]));
          v[i][j][k] = -v[i+1][j][k];
          w[i][j][k] = 0;
         }

  /* At the top-left corner-cell. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(j = NthG; j < NthG+nth; j++) u[NrG-1][j][NzG+nz] = g->u[1][NrG+g->nr[1]-g->nr_gap-1][j][NzG];

  /* At the top-interface connected to the top-section 1. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG+nz] = g->u[1][g->nr[1]-g->nr_gap+i][j][NzG];
          v[i][j][NzG+nz] = g->v[1][g->nr[1]-g->nr_gap+i][j][NzG];
          w[i][j][NzG+nz] = g->w[1][g->nr[1]-g->nr_gap+i][j][NzG];
         }

  /* At the outer side-wall. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < NthG+nth; j++)
         {
          u[NrG+nr][j][k] = (r[NrG+nr-1]/r[NrG+nr+1])*u[NrG+nr-2][j][k];
          v[NrG+nr][j][k] = -v[NrG+nr-1][j][k];
          w[NrG+nr][j][k] = -w[NrG+nr-1][j][k];
         }

  /* At the OUTLET: Orlansky BC */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic,2)
#    endif
     for(_i = 0; _i < nr; _i++) for(i = _i+NrG, j = NthG, _j = 0; _j < nth; _j++, j++)
         {
          switch(iRK)
            {
             case 0:
             u[i][j][NzG-1] = g->outBC.tnext_ug[_i][_j];
             v[i][j][NzG-1] = g->outBC.tnext_vg[_i][_j];
             w[i][j][NzG-1] = g->outBC.tnext_wg[_i][_j];
             break;

             case 1:
             u[i][j][NzG-1] = c2*(g->outBC.u[2][_i][_j][0]/6 - g->outBC.u[1][_i][_j][0] + g->outBC.u[0][_i][_j][0]/2\
                              + g->outBC.tnext_ug[_i][_j]/3) + g->outBC.u[0][_i][_j][0];

             v[i][j][NzG-1] = c2*(g->outBC.v[2][_i][_j][0]/6 - g->outBC.v[1][_i][_j][0] + g->outBC.v[0][_i][_j][0]/2\
                              + g->outBC.tnext_vg[_i][_j]/3) + g->outBC.v[0][_i][_j][0];

             w[i][j][NzG-1] = c2*(g->outBC.w[2][_i][_j][0]/6 - g->outBC.w[1][_i][_j][0] + g->outBC.w[0][_i][_j][0]/2\
                              + g->outBC.tnext_wg[_i][_j]/3) + g->outBC.w[0][_i][_j][0];
             break;

             case 2:
             u[i][j][NzG-1] = (a31+a32)*g->outBC.u[2][_i][_j][0]/6 - (a31+a32-a32*c2)*g->outBC.u[1][_i][_j][0]\
                 + (1 + 0.5*(a31+a32-4*a32*c2))*g->outBC.u[0][_i][_j][0] + (a31+a32+3*a32*c2)*g->outBC.tnext_ug[_i][_j]/3;

             v[i][j][NzG-1] = (a31+a32)*g->outBC.v[2][_i][_j][0]/6 - (a31+a32-a32*c2)*g->outBC.v[1][_i][_j][0]\
                 + (1 + 0.5*(a31+a32-4*a32*c2))*g->outBC.v[0][_i][_j][0] + (a31+a32+3*a32*c2)*g->outBC.tnext_vg[_i][_j]/3;

             w[i][j][NzG-1] = (a31+a32)*g->outBC.w[2][_i][_j][0]/6 - (a31+a32-a32*c2)*g->outBC.w[1][_i][_j][0]\
                 + (1 + 0.5*(a31+a32-4*a32*c2))*g->outBC.w[0][_i][_j][0] + (a31+a32+3*a32*c2)*g->outBC.tnext_wg[_i][_j]/3;

             break;
            }
         }

  /* It requires the outlet Orlansky-BC applied first. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(j = NthG; j < NthG+nth; j++)
         {
          u[NrG-1][j][NzG-1] = 0;
          u[NrG+nr-1][j][NzG-1] = 0;
          w[NrG-1][j][NzG-1] = -w[NrG][j][NzG-1];
          w[NrG+nr][j][NzG-1] = -w[NrG+nr-1][j][NzG-1];
         }

  /* At outlet, calculating w[i][j][NzG-2].
     It requires the outlet Orlansky-BC applied first. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic,1)
#    endif
     for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
         {
       /* Making use of the periodic-BC in the theta-direction. */
          if(j == NthG) dv = v[i][NthG][NzG]+v[i][NthG][NzG-1]-v[i][NthG+nth-1][NzG]-v[i][NthG+nth-1][NzG-1];
          else dv = v[i][j][NzG]+v[i][j][NzG-1]-v[i][j-1][NzG]-v[i][j-1][NzG-1];

          if(i == NrG) rdu = r[i+1]*(u[i][j][NzG]+u[i][j][NzG-1]); /* Because of the solid wall. */
          else if(i == NrG+nr-1) rdu = -r[i]*(u[i-1][j][NzG]+u[i-1][j][NzG-1]);  /* Because of the solid wall. */
          else rdu = r[i+1]*(u[i][j][NzG]+u[i][j][NzG-1]) - r[i]*(u[i-1][j][NzG]+u[i-1][j][NzG-1]);

          w[i][j][NzG-2] = w[i][j][NzG] + (z[NzG+1]-z[NzG-1])*(rdu/(r[i+1]-r[i]) + dv/(th[j+1]-th[j]))/(r[i]+r[i+1]);
         }

  /* For the ghost cells in theta-direction. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(dynamic,3)
#    endif
  /* It needs all the other BCs to be applied first. */
     for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
      /* Applying Periodic B.C. */
         u[i][NthG-1-j][k] = u[i][NthG+nth-1-j][k];
         u[i][NthG+nth+j][k] = u[i][NthG+j][k];

         v[i][NthG-1-j][k] = v[i][NthG+nth-1-j][k];
         v[i][NthG+nth+j][k] = v[i][NthG+j][k];

         w[i][NthG-1-j][k] = w[i][NthG+nth-1-j][k];
         w[i][NthG+nth+j][k] = w[i][NthG+j][k];
        }
   }

   return;
  }
